Network inference using asynchronously updated kinetic Ising Model 
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Network structures are reconstructed from dynamical data by respectively naive mean field (nMF) 
and Thouless- Anderson-Palmer (TAP) approximations. For TAP approximation, we use two meth- 
ods to reconstruct the network: a) iteration method; b) casting the inference formula to a set of 
cubic equations and solving it directly. We investigate inference of the asymmetric Sherrington- 
Kirkpatrick (S-K) model using asynchronous update. The solutions of the sets cubic equation 
depend of temperature T in the S-K model, and a critical temperature T c is found around 2.1. For 
T < T c , the solutions of the cubic equation sets are composed of 1 real root and two conjugate 
complex roots while for T > T c there are three real roots. The iteration method is convergent 
only if the cubic equations have three real solutions. The two methods give same results when the 
iteration method is convergent. Compared to nMF, TAP is somewhat better at low temperatures, 
but approaches the same performance as temperature increase. Both methods behave better for 
longer data length, but for improvement arises, TAP is well pronounced. 

PACS numbers: 02.50.Tt, 02.30.Mv, 89.75.Fb, 87.10.Mn 



I. INTRODUCTION 

A present challenge in biological research is how to deal 
with the data originating from the high-throughput tech- 
nologies. Information can often convincingly be struc- 
tured in the form of networks Vertices on a network 
are entities and the links with numbers or other descrip- 
tions attached to them are the interactions between the 
elements, in, e.g., the biological system [HI]- On dif- 
ferent levels of abstraction, information about the inter- 
actions between each pair of elements is hence useful to 
understand the biological system. Finding interactions 
between entities from the empirical data is an inverse 
problem called 'network reconstruction' 0, @, 0, ITTl - fl3| . 

In this work, we use an idealized system to generate 
'empirical' data with computer, and then try to recon- 
struct the network structure of the system, using this test 
data. The system is the kinetic Ising model, intended as 
a proxy for simultaneous recordings from many neurons. 
In this setting, symmetric couplings between the entities 
are not appropriate, since two neurons will typically not 
act on each other in a symmetric way fl5j . The prop- 
erties of asymmetric neural networks have been studied 
previously but not much work has been done 

in the context of network reconstruction. Here we ex- 
tend a presentl y rep orted approach using dynamic mean 
field theory [lfX l29j j from synchronously updated models 
to asynchronously updated models. The analysis closely 
parallels that of [29| , with the difference that data is con- 
tinuous in time. The similarities and differences between 
our results and [29} are commented upon in Conclusion. 

Multi-neuron firing patterns can be observed with 
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present technologies up to thousands of neurons ( record- 
ings on retina systems ). Schneidman et al. @ showed 
that the interactions between neuron pairs could be re- 
constructed using only the observed firing rates and the 
pair- wise correlations. Recently, questions have arisen 
whether the methods used in Q generalize to other data 
sets, and if the approximations involved can be improved 
or not [1, ITTl - [l3l | . There has also been significant devel- 
opment on the more theoretical side [1, [hi], H3, [H| . 

A theoretical model, which can be used to generate the 
frequencies of all possible spiking configurations is the 
well-known Ising model Q. For a system of N neurons, 
it is characterized by up to iV 2 parameters: N external 
fields, 0j, on each individual neuron, and N(N—1) 'links', 
Jij, between each pair of neurons. In the asymmetric 
Ising model, J^- is not equal to Jji. And for S-K model, 
the symmetrized and anti-symmetrized couplings J?j and 
J°j s are identically independent Gaussian distributed ran- 
dom variables. The model is entitled 'kinetic' because, 
except for the fully symmetric case, it does not corre- 
spond to an equilibrium statistical mechanics system. 

With the observed average firing rates and all pairwise 
equal-time correlations in an empirical data set, maxi- 
mum entropy models can find a probability distribution 
which maximizes the entropy of the data domain. This 
condition implies that the samples are drawn indepen- 
dently from the same distribution. The state of max- 
imum entropy given is an equilibrium state which has 
a probability distribution of Ising form Q. The quan- 
tities Jij and hi are then Lagrange multipliers to sat- 
isfy the constrains that the ensemble expectation values 
agree with sample averages in the data set. If the data is 
however generated by a dynamics, then samplings drawn 
close in time are typically dependent. This is the extra 
information which will be used here through the kinetic 
inverse Ising reconstruction scheme. For the equilibrium 
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version of the inverse Ising problem, Yasser Roudi and 
collaborators review and investigate several approxima- 
tion methods [111, EH HH with the maximum entropy 
method, arriving at the general conclusion that all of 
them are unreliable in a dynamic setting, if the systems 
are sufficiently large, and in most ranges of parameters. 
Better inference methods on dynamic data are called for. 

A standard approach to sample the equilibrium Ising 
model is Glauber dynamics, which we will describe below. 
It is however not restricted to symmetric Ising model, but 
also well-defined for models with asymmetric couplings. 
It is plausible that such a more general frame-work can 
describe the underlying system not close to equilibrium, 
and with asymmetric couplings, better. Here we are 
therefore interested in using kinetic Ising model, typi- 
cally with asymmetric couplings, to reconstruct a neural 
network dynamically. 

There are several reasons to consider asynchronous up- 
date models (Glauber dynamics) instead of synchronous 
update. The first is that asynchronous updates converge 
to a stationary state which for symmetric models in the 
Boltzmann-Gibbs equilibrium measure, while neither is 
necessarily true for synchronous updates. A second rea- 
son is that most plausible applications are naturally asyn- 
chronous. For instance, the expression of gene is not a 
synchronous process, the transcription of DNA and the 
transport of enzymes may take from milliseconds up to 
a few seconds. Another example is the refractory pe- 
riod for neuron in which the neuron cannot respond to 
input signal as it is still processing or recovering from 
the previous input signal. The period generally lasts for 
one millisecond [l9|. Besides, [2(| Hl[ show that the bi- 
ological networks do not have a completely synchronous 
update. For these reasons, we have focused on the asyn- 
chronous update Glauber dynamics. For a discussion of 
synchronous update we refer to [HI HI] ■ 

The paper is organized as follows: we describe the 
asymmetric S-K model and Glauber dynamics in Sec. 
II; the inference formula with nMF and TAP approxi- 
mation for asynchronous case is derived in Sec. Ill; the 
performances of the inference formula are given in Sec. 
IV. Finally, we summarize the work in Sec. V. 



II. ASYMMETRIC S-K MODEL AND 
GLAUBER DYNAMICS 



consists both of identically and independently Gaussian 
distributed random variables with means and variances: 



i 



N 1 



(2) 



The self-connections are avoided, i.e., the on-diagonal el- 
ements of Jf, and Jf/ equals 0. 

We now define the kinetic Ising model with asyn- 
chronous updates. Let the joint probability distribution 
of spin states in system at time t as p(si, ...,SN;t), and 
let the master equation of our model be written as 



— p(si,...,s N ;t) 

= ^ ^i{~ s t)p( s ii -Si, s N ; t) - o; 4 (si)p(s; t),(3) 



where u>i(si) is the flipping rate, i.e., the probability for 
the state of ith neuron changes from Si to — Si per unit 
time. The flipping rates are given by Glauber dynamics 
as follows: 



1 + exp[2f3s i (9 l + V JijSj 



(4) 



where, /3 is the inverse of temperature T. For conve- 



nience, define Hi = J^. JijSj 



as the effective field 



on neuron i , where 0i is the external field of spin i. If 
the couplings are symmetric (i.e., J"? = 0), then the 
steady state of the dynamics given by (3) and (4) is 
p(si, s N ) oc exp(0J2i s $i + J2ij SiSjJij). If the cou- 
plings are not symmetric, then (3) and (4) still have a 
steady state (under general condition), but this state 
does not have a simple description. 

With state for each neuron Si , we can naturally define 
the time dependent means and correlations as follows: 



m-i = (si{t)). 
Cij(t - t ) = (si(t)sj(t )) - mirrij. 



(•5) 



From equation (3) and (4), we get the equation of mo- 
tion for means and correlations as 



drrii 
~dt 



rrn + (tanh^Si-ffi^)]). 



dt 



[ Si (t) Sj (t )) = -(si(t)aj(to)) + (tanh[/3F l (i)s J (t )]).(6) 



The S-K model is a system of N spins, which models 
N neurons with binary states (s^ = 1 for firing state, 
otherwise = — 1 ). It is a fully connected model, i.e., all 
neurons in the system have interactions with each other. 
The interactions Jij between each pair of neurons have 
the following form: 



Jii J i 



kJf, 



k > 0. 



(1) 



where, k measures the asymmetric degree of these in- 
teractions, Jf, 



asymmetric matrices J^ s 



and J?f are symmetric Jh 



Jji and 



Jjf, respectively. They 



For the second equation of eq. (6), the term in the 
left hand side and the first term in the right hand side 
can be solved based on the empirical data produced by 
the Glauber dynamics. However, the calculation of the 
average value for ta.nh[f3Hi(t)sj(to)] involves all kinds of 
higher-order correlations and is therefor not easily ex- 
pressed only in terms of means and pair- wise correlations. 
In order to solve the second equation in (6), perturba- 
tively approximations for the second term of the right 
hand side are obviously needed. Here, we use the nMF 
and TAP approximations respectively to deal with this 
tanh function. 
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III. NMF APPROXIMATION AND TAP 
APPROXIMATION 

The simplest method to find out the parameters of the 
Ising model from empirical data is the mean-field theory: 



i. L = tanh/3(0, + Jij m j 



(7) 



Following recent practice, and to distinguish this first 
level of approximation from others, we will refer to it as 
naive mean- field (nMF). Let bi = 9i + . Jijrrij and 
rewrite H; as 



Hi = bj + Jij(s 



3 

3 3 



X-MX. • (8) 



Expanding the tanh function with respect to f3bi in equa- 
tion (6) 



Jij is obtained, which is formally the same as in the nMF 
approximation. 

J = TA- 1 DC~ 1 . (14) 
However, matrix A in TAP formula is different 



A; 



%(1 



1-/3 2 (1 



(15) 

Eq. (14) is a function of the couplings J, and therefor it 
is a nonlinear equation for matrix J. 

We try to solve eq. (14) for J though two approaches. 
One way is to solve it iteratively. We start from reason- 
able initial values Jfj and insert them in the right hand 
side of the formula. The resulting is the solution af- 
ter one iteration. This can be again replaced in the right 
hand side to get the second iteration results and etcetera 



-( Si (t) Sj (t )} + ( Si (t) Sj (t )) 



niimj 



13(1 - mj) (^2 J tt (M*)*«i(*o)}^ • (9) 



and denoting the time difference t — to as r, we have 
±Cii (r) + C ti (r) = /3(1 - m\) ]T J ik C kj (r). (10) 

k 

In the limit t — > 0, we obtain the equation which we need 
to infer the network couplings: 



(11) 



where D = C + C and Aij = <% (1 — mf). 

Equation (11) is a linear matrix equation with respect 
to . We can solve it directly. 

Next, we turn to derive the inference formula with TAP 
approximation. If the Onsager term, i.e., the effect of the 
mean value of neuron i on itself via its influence on an- 
other neuron j, is taken into account, the TAP equation 

is m 

mi = tanh(/3&i - m^ 2 ^ J? fc (l - ml)). (12) 

With 

Ti = h ± m^ 2 4(1 " m t) + JikSak. (13) 

and eq. (12), we expand the tanh function in eq. (6) 
with respect to 

^-m^ 2 ^;4(l-m 2 ) 

to the third order and keep the terms only up to the third 
of J. Then the corresponding TAP inference formula for 



J t + 1 = TA(J t )- 1 DC~ 



(16) 



An alternative way is solving it directly, as done for the 
synchronous update model in (29j . casting the inference 
formula to a set of cubic equations. For eq. (15), we 
denote 



F* = /3 2 (l-m 2 )£ 4(1-7^) 



(17) 



and plug it into eq. (14), and then obtain the following 
equation for Jif. 



J 



TAP 



T * Vu 



(18) 



where Vy = [DC . Inserting eq. (18) into eq. (17), 
we obtain the cubic equation for Fi as: 



Fi{l-Fif 



E,^(i 



1 — mj 



= 0. (19) 



With the obtained physical solution for Fi, we get the 
reconstructed couplings J TAP as 



J 



TAP _ J ij 



1-Fi 



(20) 



It is worth mentioning that for the cubic equation (18), 
we have three solutions with possible imaginary parts. 
Here we study the real roots of the cubic equation and 
ignore those solutions with imaginary parts. When three 
solutions are all real ones, we take the smallest one. 

We introduce A to measure the difference between the 
reconstructed network structure and the original true 
ones, i.e., A is the reconstruction error 



A 



4) 2 



E(4) 2 



where Jy represents the true network couplings and J\, 
for the reconstructed ones. 
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IV. THE PERFORMANCES OF NMF AND TAP 
APPROXIMATION 

As the starting point, we take a look at the number 
of solutions given by nMF and TAP approximation. The 
nMF gives unique solution while the iteration method of 
TAP starting from nMF provides solution when the it- 
eration is divergent and 1 solution for convergence. How- 
ever, the cubic-equation method of TAP approximation 
always contains at least one solution. Denote the con- 
stant term of eq. (19) as x, 



x is temperature dependent and negative as < mf < 1. 
The cubic equation (19) has 3 real roots when — ^ < 
x < 0. We only consider the smallest one and indeed it 
provides the most accurate Jy's (data are not shown). 
With x < —-§f, eq- (19) has only one real root and 
other two complex solutions with imaginary part which 
are discarded as they have no physical meaning. In Fig. 
1, we give the fraction of cubic-equation set (19) (as i — 
1,2,...,N, where N is the system size) which contains 
three real solutions. When the set of cubic-equation at 
given T contains N real and 2 * N complex solutions, 
we say the fraction of three real roots equals at this 
temperature point. As shown in Fig. 1, a transition 
seems to occur around T c = 2.1. For large system size 
and T < 2.1, the solutions for eq. (19) has only one real 
root while for T > 2.1 3 real ones. We plot this figure 
for data length L — N * 10 6 , so smaller N means shorter 
data length, that explains why the curve of N = 20 is 
not quite smooth. 




FIG. 1: (Color online) The fraction of 3 real roots for the 
cubic equation set eq. (19). A transition seems to occur 
around T c = 2.1. Here, we find larger N, the transition curve 
sharper. The parameter values: 9 — 0.5, k — 1, L — 20 * 10 6 . 

For the simulation of the iteration method of TAP ap- 
proximation, we take the reconstructed Jff F by nMF as 



the initial input J?-, and follow eq. (16) to get J\p Jfj — 
iteratively. If the average value of S(t) = | Jy — j\j \ 
less than the threshold value 10 -5 , then, we consider the 
iteration is convergent and stop iteration. An interest- 
ing phenomena of the iteration method is it is divergent 
when the solutions of cubic-equation set contain complex 
roots while convergent when they contain only real roots. 
Here, we mention three possible causes for the divergence. 
One originates from the frozen states of spin-glass where 
mf = 1 and neither nMF nor TAP can work. A second 
possible cause: there exists a single fixed point of the so- 
lution but the initial Jy 's are drawn as Jf^ F , which may 
a little bit far away from the true solutions for Jy's at 
low T, and the iteration can not reach to the fixed point. 
The last possible cause may come from the fixed point 
which is unstable. Here, the given results are for 8 = 
and k = 1, there is no frozen states for the given temper- 
atures. Then, the divergence may arise by the second or 
third possible reason. 




FIG. 2: (Color online) The reconstruction error A with tem- 
perature T for both nMF and TAP approximation. The other 
parameter values: TV = 20, temperature L — 20 * 10 10 , exter- 
nal field 8 = 0, asymmetric degree k = 1. Notations: black 
square for nMF red circle for cubic equation method for TAP, 
blue triangle for iteration method for TAP. Each data point 
is averaged on 10 realizations. 

We next turn to investigate the influence of T on the 
reconstruction errors A in the case of zero external field 
(6> = 0) aS-K model (k = 1). We plot A with T for nMF 
and TAP in Fig. 2. For TAP approximation, when iter- 
ation method is convergent, it produces the same results 
(blue triangle) as the cubic-equation method (red circle) . 
Both approximations work better with temperature T 
increasing but approach to the same behavior when T 
goes higher. It is because for eq. (15), the Onsager term 
will approach if T goes high enough, i.e., there will be 
no difference between nMF and TAP approximation. As 
shown in Fig. 2, TAP always works better than nMF 
before they approach to the same results. But there is 
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an noticeable area in which the curve by the cubic equa- 
tion method of TAP pointed to with letter 'A' is not as 
smooth as that of nMF. The reason is this temperature 
interval is located in the critical area where the solutions 
of the cubic-equation set eq. (19) are coexistence of two 
states: some spins have 3 reals roots and the others have 
only 1 real root. We tested also for systems with different 
size and found that larger system size give more clear in- 
flexions and closer to the critical temperature T c , around 
2.1. Such results are consistent with the results shown in 
Fig. 1. 

Fig. 3 illustrates the reconstructions errors for every 
Jij 's with scatter plots of the inferred Jij 's by nMF and 
TAP approximation against J*™ e 's. The left plot is for 
the data length L = N* 10 5 and L = iV* 10 7 for the right 
one. Here, the system size N = 20 and the temperature 
T = 3.7 for this plot where the iteration method of TAP 
is convergent. The scatter plot shows that both nMF and 
TAP perform better for larger L. As shown in both left 
and right hand side of Fig. 3, the data points for J™ p 's 
inferred by cubic-equation method are almost covered by 
that for Jjf P, s inferred by iteration method, especially 
for L = N * 10 7 . 




-0.6 -0.3 0.0 0.3 0.6 -0.6 -0.3 0.0 0.3 0.6 
.true 



FIG. 3: (Color online) The scatter plot for the reconstructed 
couplings versus the true ones. The parameter values: N = 
20, temperature T = 3.7, external field 9 = 0, asymmetric de- 
gree k — 1. Notations: black square for inferred couplings 
using nMF versus J*™ 6 , blue circle for iteration equation 
method of TAP versus Jij Ue , red triangle for cubic method 
of TAP versus J*j' uo . 

From the right plot in Fig. 3, it is difficult to say 
which approximation is better as the reconstruction error 
is quite small especially with longer data length. Thus, 
we move next to see how the data length L works on the 
reconstruction error A in the case of zero external field 
aSK model. With the asynchronously updating Glauber 
dynamics, longer data length L (L = N * L' , where V is 
the data length in the corresponding synchronous update 
case, N is the system size) is needed to obtain compa- 
rable results with that in synchronous case [29| and say 
something about our system. In Fig. 4, A versus L for 
both nMF and TAP are plotted for a given temperature 
T = 8, where the iteration method of TAP is conver- 
gent. They both reconstruct better with increasing L, 
i.e., A decreases as L increases. For short data length 




"| Q - 1 1 1 1 1 1 i 

10 5 10 s 10 7 10 8 10 9 10 10 10 11 
L 

FIG. 4: (Color online) Reconstruction error A versus the data 
length L for nMF and TAP approximation. The other param- 
eter values: N = 20, temperature T = 8, external field 9 = 0, 
asymmetric degree k = 1. Notations: black square for nMF, 
blue circle for iteration equation method of TAP, red triangle 
for cubic method of TAP. 

L < N * 10 7 , nMF and TAP produce almost the same 
reconstruction error. However, TAP works better than 
nMF when L > N * 10 s . The A for TAP is one order 
smaller that that for nMF when L > N * 10 9 . Here, 
again, we find the data points for cubic-equation method 
of TAP are covered by those for iteration method of TAP. 

The above results are general to different system size 
N. The performances for nMF and TAP are also com- 
pared with non-zero external field 9^0. We find there 
exists a frozen state in the low testing temperature where 
neither nMF nor TAP can work there. 



V. CONCLUSION 

We studied the network inference using asynchronously 
updated kinetic Ising model. Two approximations, nMF 
and TAP, are introduced to infer the connections and 
connection strengths in the network. We have found the 
transition of the solutions' type for the cubic equation 
method of TAP with critical temperature T c w 2.1. We 
have implemented the TAP approximation as two differ- 
ent schemes, the cubic scheme, and the iteration scheme. 
For large system, the T c seems to be the starting tem- 
perature point for TAP iteration method to converge. 

Comparing our work with [29| in which the syn- 
chronously updated Glauber dynamics is used, we find 
two similarities. The first one is both approximations 
reconstruct better with increasing temperature or longer 
data length. The other one is TAP works better than 
nMF especially with long data length at given tempera- 
tures. There are also differences. For instance, the im- 
provement by TAP approximation in asynchronous case 
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is not as much as that in synchronous case. Besides, 
in order to get the comparable results with synchronous 
case, the data length for asynchronous case should be at 
least N times longer than that for synchronous case. 

This work is able to extend to deal with the biological 
data from experiments, especially for data produced in 
continuous time which correspond to the asynchronous 
updates. Given the large amount of data needed to see a 
difference, we believe that in most application scenarios, 
network inference using asynchronously updated kinetic 
Ising models should work well enough using naive mean- 
field (nMF) reconstruction, and the further step to TAP 
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